2008-01-0777 


Experimental  Validation  of  a  Coupled  Fluid-Multibody 

Dynamics  Model  for  Tanker  Trucks 

Tamer  M.  Wasfy  James  O’Kins;  Scott  Smith 

Advanced  Science  and  Automation  Corp.  U.S.  Army  Tank  Automotive  Research  Development  and 

Indianapolis,  IN  Engineering  Center,  Warren,  Ml 

Copyright  ©  2008  SAE  International 


ABSTRACT 

A  time-accurate  finite  element  model  for  predicting  the 
coupled  dynamic  response  of  tanker  trucks  and  their 
liquid  payloads  is  presented  along  with  an  experimental 
validation  of  the  model.  The  tanker  truck  components 
are  modeled  using  rigid  bodies,  flexible  bodies,  joints 
and  actuators.  The  model  is  validated  using  a  full-scale 
army  heavy  class  tactical  trailer  carrying  a  water  tank. 
The  trailer  is  placed  on  an  n-post  motion  base  simulator 
which  was  used  to  perform  harmonic/ramp  pitch,  roll  and 
stir  excitations  of  the  trailer  in  order  to  simulate  typical 
road  maneuvers.  Experiments  were  carried  out  with  an 
empty  tank  and  a  65%-filled  tank.  The  time-histories  of 
the  tires  and  suspension  system  deflections  are 
measured  for  the  various  input  motion  excitations.  The 
experiment’s  measurements  are  compared  with  the 
results  predicted  using  the  computational  model.  The 
comparison  shows  that  the  model  can  predict  with 
reasonably  good  accuracy  the  test  tanker-trailer’s 
dynamic  response. 

1.  INTRODUCTION 

A  time-accurate  finite  element  model  for  predicting  the 
coupled  dynamic  response  of  tanker  trucks  and  their 
liquid  payloads  is  presented  along  with  an  experimental 
validation  of  the  model.  Tanker  trucks  are  used  to 
transport  liquids  (such  as  fuels  and  water)  from 
production  facilities  to  distribution  outlets.  They  provide 
the  most  flexible  and  readily  deployable  form  of  liquid 
transportation.  Tanker  trucks  usually  consist  of  a  truck 
that  may  be  pulling  one  or  more  trailers  with  one  or  more 
tanks  carried  on  the  truck  and/or  the  trailer(s).  Figure  1 
shows  a  typical  military  tanker  truck  system  consisting  of 
a  truck  pulling  a  trailer  that  is  carrying  a  potable  water 
tank. 

Many  tanker  truck  accidents  occur  due  to  sloshing  of  the 
liquid  payload  when  the  tank  is  nearly  half-full.  The 
sloshing  can  shift  the  center  of  mass  of  the  tank  and  can 
cause  high  time-varying  inertia  forces.  This  can  cause 
the  tanker  truck  to  overturn,  lose  traction  on  one  or  more 


tires,  or  be  hard  to  brake  or  steer.  All  this  can  result  in 
tanker-truck  accidents  which  are  costly  and  dangerous 
especially  if  the  tank  contains  hazardous  liquids  (such  as 
fuels,  petrochemicals,  or  toxic  liquids.)  The  estimated 
yearly  cost  of  truck  accidents  in  the  US  is  about  $1 .2 
billion  for  trucks  carrying  hazardous  material  and  $43 
billion  for  trucks  carrying  non-hazardous  material  [1]. 
The  majority  of  those  accidents  (more  that  75%)  occur 
en-route.  Accidents  involving  fire  trucks  account  for 
approximately  20%  of  U.S.  fire  fighter  deaths  each  year. 
Crashes  involving  fire  trucks  are  the  most  prevalent  of 
these  types  of  accidents  with  87%  occurring  when  the 
fire  truck  went  out  of  control  and  rolled  over  or  left  the 
road  (no  collision)  [2].  The  cost  of  physically  building  and 
testing  tanker  trucks  can  be  very  high.  Several 
configurations  with  various  shapes  and  locations  for  the 
baffles  within  the  tanks,  different  tank  geometries,  and 
suspension  system  characteristics  need  to  be  tested. 
Considerable  savings  and  better  performance  can  be 
achieved  if  most  of  the  testing  is  done  numerically.  The 
model  presented  and  validated  in  this  paper  can  help 
increase  the  stability,  safety  margins,  maximum 
allowable  speed,  and/or  liquid  carrying  capacity  of  tanker 
trucks  and  can  help  to  reduce  the  annual  costs  of  tanker 
truck  accidents.  It  can  also  be  used  to  investigate  the 
cause  of  tanker  truck  accidents  and  to  come  up  with  the 
most  effective  design  improvements  that  can  increase 
the  safety  of  tanker  trucks. 


Figure  1  An  FMTV  (Family  of  Medium  Tactical  Vehicles)  military  tanker 
truck  system  consisting  of  a  truck  pulling  a  trailer  that  is  carrying  a 
water  tank. 


The  model  can  also  be  easily  extended  to  other 
applications  involving  flexible  multibody  systems 
carrying  liquid  filled  tanks  such  as:  tanker  railroad  cars, 
tanker  ships,  passenger  cars  (fuel-tanks),  airplanes 
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(fuel-tanks  and  payload  tanks),  space  launch  vehicles 
(liquid  fuel  tanks),  and  satellites/space  structures  (liquid 
fuel  tanks). 

The  experimental  validation  of  the  numerical  model  was 
carried  out  using  a  full-scale  PLS  (Palletized  Load 
System)  tanker-trailer  that  was  instrumented  and 
mounted  on  an  n-post  motion  simulator  (Figure  2).  The 
n-post  motion  simulator  consists  of  hydraulic  linear 
actuators  under  each  tire  that  can  be  independently 
controlled  to  simulate  typical  road  maneuvers  that  the 
vehicle  will  be  subjected  to  such  as  lane-change  and 
going  over  symmetric  or  asymmetric  bumps/potholes. 
Experiments  were  carried  out  with  an  empty  tank  and  a 
65%  filled  tank.  The  time-histories  of  the  tires  and 
suspension  system  deflections  are  measured  for  various 
input  motion  excitations.  The  measurements  are 
compared  with  the  results  predicted  using  the 
computational  model. 


Figure  2  Instrumented  PLS  (Palletized  Load  System)  tanker-trailer 
mounted  on  an  n-post  motion  simulator  that  was  used  to  validate  the 
present  finite  element  model. 


In  order  to  accurately  predict  the  dynamic  response  of 

tanker  trucks,  the  model  must  accurately  account  for  the 

following  effects: 

•  Incompressible  liquid  flow  in  a  moving/deforming 
container  (governed  by  the  incompressible  Navier- 
Stokes  equations). 

•  Modeling  of  the  liquid’s  free-surface. 

•  Modeling  of  fluid  viscosity  and  turbulence. 

•  Coupling  between  the  solid  and  the  fluid  at  the  fluid- 
structure  interface. 

•  Large  rotation  of  the  solid  bodies. 

•  Deformation  of  the  solid  flexible  bodies.  Flexible 
bodies  can  be  modeled  as  springs,  beam,  shells  or 
general  solids.  For  example  a  general  solid  model 
can  be  used  for  leaf-springs  while  a  linear-spring 


model  can  be  used  to  model  struts  and  suspension 
springs. 

•  Joint  kinematic  constraints  for  spherical,  revolute, 
cylindrical  and  prismatic  joints  including  joint  friction 
and  clearances. 

•  Tire  load-deflection  and  damping  characteristics. 

•  Frictional  contact  between  the  tire  and  the  terrain. 

•  Actuators  and  control  laws. 

•  Motion  control  components  including  transmission 
components,  clutches,  and  brakes  (All  those 
components  involve  friction). 

In  the  present  model,  the  fluid  governing  equations  of 
motion  are  the  incompressible  Arbitrary  Lagrangian- 
Eulerian  Navier-Stokes  equations  along  with  a  large- 
eddy  simulation  (LES)  turbulence  model.  Incompressible 
fluid  flow  can  be  modeled  using  finite  volume,  finite 
element  or  particle  methods.  In  addition,  a  technique  for 
modeling  the  free-liquid  surface  must  also  be  used.  A 
literature  review  of  the  techniques  for  modeling  fluid  flow 
with  a  free-surface  was  presented  in  [3],  Those 
techniques  are: 

(a)  Volume-of-fluid  (VOF)  method  [4-9],  Each  element 
has  a  VOF  value  between  0  (for  empty  elements) 
and  1  (for  elements  completely  filled  with  fluid).  The 
free  surface  is  reconstructed  for  each  element  using 
piecewise-linear  planar  segments  that  are  calculated 
from  the  VOF  value  of  the  element  along  with  the 
VOF  values  of  neighboring  elements  (which  are 
used  to  determine  the  normal  to  the  planar  surface). 

(b)  Level-set  method.  This  method  uses  a  smooth 
scalar  function  defined  at  every  node  in  the  fluid 
domain  which  specifies  the  signed  smallest 
Eucledian  distance  between  the  node  and  the 
interface  [10].  The  evolution  of  the  scalar  function  is 
governed  by  a  convection  transport  equation  where 
the  interface  is  moved  with  the  fluid  velocity. 

(c)  Arbitrary  Lagrangian-Eulrian  (ALE)  method.  Using 
this  method,  the  fluid  mesh  deforms  and  moves 
along  with  the  fluid’s  free-surface  [11-16],  A 
disadvantage  of  this  method  is  that  it  does  not  allow 
large  surface  deformation  including  surface  break¬ 
up  and  merging  unless  the  fluid  domain  is  re¬ 
meshed  [14,  15]. 

(d)  Lagrangian  particle  methods.  Lagrangian  particles 
which  represent  small  packets  of  fluid  are  used  to 
model  the  fluid  flow.  A  contact  model  between  the 
fluid  particles  is  used  to  model  the  fluid 
compressibility  and  viscous  effects.  A  special  type  of 
this  class  of  methods  which  has  been  successfully 
applied  to  free-surface  flows  and  fluid-structure 
interaction  problems  is  the  particle  finite  element 
method  (PFEM)  [17]  in  which  the  particles  are  used 
to  generate  a  polyhedral  finite  element  mesh  every 
time  step  using  an  extended  Delaunay  tesselation. 
The  solution  of  the  incompressible  Navier-Stokes 
equations  is  then  carried  out  using  that  mesh. 

Techniques  to  handle  fluid  flow  in  a  moving/deforming 
container  were  also  reviewed  in  [3],  They  include: 


1.  Fixed  Cartesian  fluid  mesh  with  cut-cell  boundary 
condition.  The  fluid  domain  is  a  Cartesian  mesh  [4, 
5].  The  container  moves  inside  this  mesh.  A  cut-cell 
technique  is  used  to  find  where  the  boundary  of  the 
container  intersects  with  the  fluid  cells. 

2.  Moving  ALE  mesh.  The  fluid  is  modeled  using  a  fluid 
mesh  that  moves  and  deforms  with  the  tank  [3,  9, 
11-13]. 

3.  Fixed-fluid  mesh  with  the  Navier-Stokes  equations 
written  in  a  reference  frame  fixed  to  the  tank  [11, 
18]. 

4.  Particle  methods  [17].  Normal  contact  constraint 
between  the  particles  and  the  tank  wall  is  used  to 
model  the  wall  impenetrability  constraint.  Tangential 
friction  between  the  particles  and  the  tank  wall  is 
used  to  model  wall  adhesion  and  viscous  effects. 

In  order  to  model  fluid  flow  with  a  free-surface  in  a 
moving  deforming  container,  the  above  methods  must 
be  combined.  Table  1  shows  the  references  where  the 
various  combinations  of  the  above  techniques  were 
used. 

Table  1  References  for  the  various  combinations  of  techniques  for 
modeling  a  free-surface  flow  in  a  moving  deforming  container. 


In  the  present  paper,  an  acceptor-donor  VOF  based 
algorithm  is  used  to  model  the  free-surface  and  the  ALE 
method  is  used  for  modeling  fluid  flow  in  a 
moving/deforming  container.  Using  the  VOF  technique 
for  modeling  the  free-surface  has  the  following 
advantages  over  other  methods: 

•  The  VOF  technique  can  conserve  mass  and  energy 
better  than  the  level-set  method.  This  is  due  to  the 
fact  that  the  VOF  technique  explicitly  balances  local 
mass  transfers  between  elements.  On  the  other 
hand,  the  level-set  method  relies  on  the  fluid  flow 
solution  to  conserve  mass  and  uses  interpolated 
velocities  which  can  lead  to  accumulation  of  mass 
loss  or  gain  errors. 

•  Unlike  the  ALE  technique,  the  VOF  technique  can 
handle  fluid  breakup  and  merging  without  the  need 
for  re-meshing.  If  re-meshing  is  used  frequently,  it 
can  degrade  the  solution  accuracy  due  to  re¬ 
interpolation  of  the  solution  field  onto  the  new  mesh. 

•  Particle  methods  require  a  large  number  of  particles 
to  accurately  model  the  free  surface.  The  PFEM 
requires  less  fluid  particles,  however,  the  tessellation 
step  is  computationally  intensive. 


Using  the  ALE  method  for  modeling  the 
moving/deforming  container  has  the  following 
advantages  over  other  methods: 

•  In  the  fixed  Cartesian  fluid  mesh  method,  since  the 
solid  container  cuts  the  mesh  cells  at  arbitrary 
surfaces,  the  fluid-solid  impenetrability  and  no-slip 
boundary  conditions  are  satisfied  only  in  a  time 
average  sense.  Also,  the  method  has  stability  and 
accuracy  problems  when  the  cut-cell  elements  at  the 
solid-fluid  interface  become  small.  In  addition,  for 
tanker  truck  applications,  since  the  tank  can 
undergo  very  large  translation  as  the  truck  is  moving 
on  the  road,  the  Cartesian  mesh  must  be  moved 
with  the  tank.  This  means  the  fluid  mesh  is  no  longer 
fixed,  and  that  an  ALE  formulation  must  be  used  to 
move  the  mesh.  This  reduces  the  simplicity  of  the 
method. 

•  The  main  disadvantage  of  using  a  fixed-fluid  mesh 
with  the  Navier-Stokes  equations  written  in  a 
reference  frame  fixed  to  the  tank  is  that  since  the 
tank  frame  is  a  non-inertial  frame  (accelerating 
frame),  writing  the  equations  of  motion  with  respect 
to  that  frame  results  in  a  complex  inertia  operator 
which  involves  centrifugal  and  Coriolis  acceleration 
terms.  Also,  this  method  cannot  -  by  itself  -  deal  with 
a  deforming  container. 

A  review  of  multibody  dynamics  modeling  techniques 
including  deformation  reference  frames,  treatment  of 
large  rotations,  discretization  techniques,  finite  elements, 
constraint  and  contact  modeling,  and  solution 
techniques  is  presented  in  [19].  In  the  present  paper,  a 
flexible  multibody  dynamics  code  with  the  following 
characteristics  is  used: 

•  The  solid,  fluid  and  fluid-solid  interface  governing 
equations  of  motion  are  solved  along  with 
joint/constraint  equations  using  a  time-accurate 
explicit  solution  procedure  that  can  maintain  time- 
accuracy  over  long  simulation  times.  The  explicit 
solution  procedure  was  presented  in  Ref.  [3]  and,  for 
the  sake  of  completeness,  is  outlined  in  the 
Appendix. 

•  Total  Lagrangian,  total  displacement  equations  of 
motion  formulation  with  the  degrees  of  freedom 
referred  to  a  global  inertial  reference  frame  [15-18], 

•  Flexible  bodies  can  be  modeled  using  a  library  of 
spring,  truss,  beam,  and  solid  nonlinear  finite 
elements  with  Cartesian  coordinate  degrees  of 
freedom.  This  allows  arbitrarily  large  element 
rotations.  The  elements  library  includes: 

o  Struts  and  suspension  springs,  that  are 
modeled  using  a  linear  spring  element,  that 
can  include  a  prescribed  stiffness-deflection 
relation,  a  damping-deflection  rate  relation 
and  a  friction  coefficient-deflection  rate 
relation. 

o  Torsional-spring  type  3-node  beam 
elements  [20,  21]. 

o  Natural-modes  eight-node  brick  elements 
[22,  23],  Those  elements  can  also  be  used 


to  model  shells  and  beams.  One  element 
through  the  thickness  is  sufficient  to 
accurately  model  the  membrane,  shear,  and 
bending  characteristics.  They  do  not  exhibit 
locking  or  spurious  modes  (widely  used 
techniques  to  alleviate  locking  such  as 
hourglass  control  lead  to  elements  that  do 
not  maintain  solution  accuracy  over  very 
long  solution  times).  They  are 
computationally  efficient.  Assumed  strain 
elements  of  comparable  accuracy  are  more 
computationally  expensive.  Any  material  law 
can  be  used  with  those  elements  including: 
linear  elastic,  hyper-elastic,  and  non-linear 
laws.  Leaf-springs  are  modeled  using  the 
natural-modes  brick  elements. 

•  Connection  points  on  the  rigid  and  flexible  bodies 
are  used  to  define  joints  between  the  bodies 
including  spherical,  revolute,  cylindrical,  prismatic 
and  rigid  joints.  A  penalty  model  is  used  to  impose 
the  joint/contact  constraints  [24],  An  asperity-based 
friction  model  is  used  to  model  joint/contact  friction 
[25], 

•  Rigid  body  rotational  equations  of  motion  are  written 
in  a  body  (material)  frame,  with  the  resulting 
incremental  rotations  added  to  the  total  body  rotation 
matrix  [24].  The  rigid  bodies  rotational  equations  of 
motion  are  written  in  a  body-fixed  frame  with  the 
total  rigid-body  rotation  matrix  updated  each  time 
step  using  incremental  rotations. 

•  Normal  contact  is  modeled  using  a  penalty 
formulation  [26,  27].  Frictional  contact  is  modeled 
using  an  accurate  and  efficient  asperity-based 
friction  model  [25].  Tires  are  modeled  using  a  one 
node  tire  model  that  accounts  for  distributed  tire 
normal  contact  using  the  penalty  formulation  and 
friction  using  the  asperity-based  friction  model. 

•  Special  elements  are  used  to  model  wheels/pulleys 
[26,  27],  sprockets  [28],  and  clutches  [29], 

•  The  code  uses  a  general  contact  search  algorithm 
that  finds  the  contact  penetration  between  finite 
elements  and  other  elements  as  well  as  general 
triangle  and  quadrilateral  surfaces. 

Two-way  coupling  between  the  multibody  system 
(vehicle)  motion  and  the  fluid  is  achieved  by  satisfying 
the  following  conditions  at  the  solid-fluid  interface: 

•  The  fluid  velocity  normal  to  the  solid’s  surface  must 
be  equal  to  the  normal  solid  velocity. 

•  The  fluid  velocity  tangent  to  the  solid  surface  can 
range  from  being  equal  to  the  tangential  velocity  of 
the  solid  surface  (no  slip  condition)  to  being  free. 

•  No  additional  energy  or  momentum  to  the  system 
should  be  introduced  at  the  interface. 

The  above  boundary  conditions  are  satisfied  by  solving 
Newton’s  equations  of  motion  at  the  fluid-solid  interface 
nodes  for  a  common  normal  acceleration  for  the  fluid 


and  the  solid  at  the  interface.  The  tangential  fluid  and 
solid  accelerations  can  range  from  being  the  same  (no¬ 
slip  condition)  to  being  completely  decoupled. 

In  the  present  paper,  a  single  computational  code  which 
uses  a  time-accurate  explicit  solution  procedure  is  used 
to  solve  both  the  solid  and  fluid  equations  of  motion. 
Many  commercial  software  and  studies  on  modeling 
liquid  sloshing  coupled  with  solid  body  motion  use  two 
codes  which  pass  the  interface  forces  and  motion  back 
and  forth  and  iterate  on  the  two  codes  until  equilibrium  is 
achieved  [e.g.  6,  18].  This  approach  adds  extra 
computational  burden  and,  in  general,  does  not  achieve 
the  same  accuracy  as  the  single  integrated  code 
solution  due  to  the  difficulty  in  achieving  an  equilibrium 
solution  between  two  disjointed  codes. 

The  rest  of  this  paper  is  organized  as  follows.  In  Section 
2  the  equations  of  motion  for  the  solid,  fluid  and  fluid- 
structure  interface  are  presented.  In  Section  3  the 
frictional  contact  model  is  presented  along  with  a 
description  of  the  1-node  distributed  contact  tire  model. 
In  Section  4  the  VOF  free-surface  model  is  presented.  In 
Section  5  the  experimental  setup  along  with  the 
validation  study  are  presented.  Finally,  in  Section  6 
concluding  remarks  are  offered. 

2.  EQUATIONS  OF  MOTION 

In  the  subsequent  equations  the  following  conventions 
will  be  used: 

•  The  indicial  notation  is  used. 

•  The  Einstein  summation  convention  is  used  for 
repeated  subscript  indices  unless  otherwise  noted. 

•  Upper  case  subscript  indices  denote  node  numbers. 

•  Lower  case  subscript  indices  denote  vector 
component  number. 

•  The  superscript  denotes  time. 

•  A  superposed  dot  denotes  a  time  derivative. 

2.1  SOLID  EQUATIONS  OF  MOTION 

The  solid  translational  equations  of  motion  are  written 
with  respect  to  the  global  inertial  reference  frame  and 
are  obtained  by  assembling  the  element  equations.  The 
finite  elements  used  here  use  only  translational  DOFs 
with  no  rotational  DOFs.  The  translational  equations  of 
motion  include  the  rigid-body  translational  DOFs.  They 
can  be  written  as: 

MKXtKi=FstKi+FatKi  (1) 

where  t  is  the  running  time,  K  is  the  global  node  number 
(no  summation  over  K;  K=1->A/  where  N  is  the  total 
number  of  nodes),  /  is  the  coordinate  number  (/=1 ,2,3),  a 
superposed  dot  indicates  a  time  derivative,  MK  is  the 
lumped  mass  of  node  K,  x  is  the  vector  of  nodal 
Cartesian  coordinates  with  respect  to  the  global  inertial 
reference  frame,  and  x  is  the  vector  of  nodal 
accelerations  with  respect  to  the  global  inertial  reference 
frame,  Fs  is  the  vector  of  internal  structural  forces,  and 
Fa  is  the  vector  of  externally  applied  forces,  which 
include  surface  forces  and  body  forces. 


For  each  rigid  body,  a  body-fixed  material  frame  is 
defined.  The  rigid  body  is  represented  by  one  node 
located  at  the  body’s  center  of  mass,  which  is  also  the 
origin  of  this  frame.  The  mass  of  the  body  is 
concentrated  at  the  node  and  the  inertia  of  the  body, 
given  by  the  inertia  tensor  ltj ,  is  defined  with  respect  to 
the  body  frame.  The  orientation  of  the  body-frame  is 
given  by  R‘°  which  is  the  rotation  matrix  relative  to  the 

global  inertial  frame  at  time  The  rotational  equations 
of  motions  are  written  for  each  rigid  body  with  respect  to 
its  body-fixed  material  frames  as: 

hAj  =  T'n  +  Ta>Ki  ~  fe  ><  VkAj ))Ki  (2) 

where  IK  is  the  inertia  tensor  of  rigid  body  K,  q  and  q 

are  the  angular  acceleration  and  velocity  vectors’ 
components  for  rigid  body  K  relative  to  it’s  material 
frame  in  direction  j,  TsKi  is  the  component  of  the  vector  of 
internal  torque  at  node  K  in  direction  i,  and  TaKi  is  the 
component  of  the  vector  of  applied  torque  in  direction  i. 
The  summation  convention  is  used  only  for  the  lower 
case  indices  i  and  j. 

The  trapezoidal  rule  is  used  as  the  time  integration 
formula  for  solving  equation  (1)  for  the  global  nodal 
positions  x: 

x'kj  =  x‘k +  0-5  A  t  (x‘Kj  +  x‘kja‘)  (3a) 

xkj  =  x ‘k/“ '  +  0.5  Ar  (x'Kj  +  x‘KjAi  )  (3b) 

where  At  is  the  time  step.  The  trapezoidal  rule  is  also 
used  as  the  time  integration  formula  for  the  nodal 
rotation  increments: 

9^=9^' +0.5  Atid^+e1-;1)  (4a) 

A&Kj  =  0.5  A?  (0'Kj  +  (4b) 

where  A()K]  are  the  incremental  rotation  angles  around 
the  three  body  axes  for  body  K.  The  rotation  matrix  of 
body  K (Rk)  is  then  evaluated  using: 

R'k  =  R(A0'Ki)  (5) 

where  R(AffKj)  is  the  rotation  matrix  corresponding  to 

the  incremental  rotation  angles  from  Equation  (4b). 

The  explicit  solution  procedure  used  for  solving 
equations  (1-5)  along  with  the  constraint  equations  is 
presented  in  Section  6.  The  constraint  equations  are 
generally  algebraic  equations,  which  describe  the 
position  or  velocity  of  some  of  the  nodes.  They  include: 

•  Prescribed  motion  constraints: 


The  dynamic  response  of  the  fluid  is  described  by  the 
ALE  version  of  the  incompressible  Navier-Stokes 
equations.  Namely,  the  equations  of  conservation  of 
momentum  and  mass  for  a  moving  deforming  control 
volume  along  with  a  large-eddy  simulation  (LES) 
transport  equation.  Those  are: 

v  «  V  «  V  VKj  V 

r±dv+rf(P^)dv  =  0  (10) 

J  /9i  i  Pr 


p  =  p0  +  rP 


u,  =  u,  -  V, 


*,j  =  *DkkS,j  +  2(A  +  P-)D„ 

Dtj  =  O^idUf/dCj  +  diij /dCj) 


J  Fit 


(piifK) 


r  v\pu,i 

l  Sxi 

if  2/*  A 


£  =  1.0  Kiu 
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(a  +  a)- 
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where  V  is  the  element  volume,  t  is  the  running  time,  p  is 
the  density  of  the  fluid,  u  is  the  fluid  velocity  vector 
relative  to  the  global  reference  frame,  h  is  the  fluid 
velocity  vector  relative  to  the  moving  fluid  mesh,  v  is  the 
velocity  of  the  fluid  mesh,  P  is  the  relative  pressure,  x  is 
the  position  vector,  r  is  the  deviatoric  stress  tensor,  D  is 
the  rate  of  deformation  tensor,  /  is  the  body  force 


vector,  r  is  the  artificial  compressibility  parameter,  p0  is 
the  nominal  fluid  density,  p  is  the  fluid  viscosity,  pi  is  an 
additional  turbulence  viscosity  calculated  using  Equation 
15b),  K  is  the  eddy  kinetic  energy,  and  s  is  the  sub-grid 
scale  eddy  kinetic  energy  dissipation  term. 
Incompressible  flow  is  modeled  using  the  artificial 
compressibility  technique  [30],  A  finite  element 
formulation  is  used  to  derive  the  element’s  semi-discrete 
equations  of  motion  from  the  governing  equations  (9- 
13).  8-node  hexahedral  elements  are  used  with  tri-linear 
equal-order  velocity  and  pressure  interpolation.  A 
pressure  averaging  algorithm  [31]  is  used  to  eliminate 
pressure  checker-boarding  (due  to  the  use  of  an  equal 
order  interpolation  for  pressure  and  velocity).  The 
element  equations  are  assembled  into  the  global  semi¬ 
discrete  equations  of  motion: 


/({*}, 0  =  0 

(6) 

l  i  .  Vi  _  Ffm 

(16) 

•  Joint  constraints: 

VfN  P'm  =  Qfh 

(17) 

/(M)  =  0 

(7) 

II 

-  ^ 

(18) 

•  Contact/impact  constraints: 

/({*})  >0  (8) 

The  penalty  technique  is  used  for  imposing  the 
constraints  in  which  a  normal  reaction  force  is  generated 
when  a  node  penetrates  into  a  contact  body.  The 
magnitude  of  the  force  is  proportional  to  the  penetration 
distance  [24-26]. 


where  MfN  is  the  lumped  fluid  mass  of  node  N,  uNi  is 
component  /  of  the  fluid  acceleration  at  node  N,  K  is 

component  /  of  the  fluid  forces  at  node  N,  VfN  is  the 
lumped  fluid  volume  at  node  N,  p‘N  is  the  fluid  pressure 

rate  at  node  N,  Qf‘N  is  the  fluid  pressure  flux  at  node  N, 
K'n  is  the  eddy  kinetic  energy  rate  at  node  N ,  and  S/N  is 


the  eddy  kinetic  energy  flux  at  node  N.  Those  equations 
are  integrated  using  the  trapezoidal  rule  along  with  an 


2.2  FLUID  EQUATIONS  OF  MOTION 


explicit  solution  procedure  to  yield  the  nodal  fluid  velocity 
and  pressure: 


ukj  ~  ukjA>  +  0-5  A t  (u‘Kj  +  u‘Kja>) 

(19) 

P^=P^+0.5At(P‘J+P‘jA‘) 

(20) 

K‘Kj=K‘ka'  +0.5  At  (K'Kj+K'kAi) 

(21) 

2.3  FLUID-STRUCTURE  INTERFACE  EQUATIONS  OF 
MOTION 

Newton’s  equations  of  motion  are  used  to  find  a 
common  normal  acceleration  for  the  fluid  and  the  solid  at 
the  interface.  This  is  done  for  each  node  at  the  fluid- 
structure  boundary  as  follows: 

(ms  +  m/)un  =  y  Fluid  Forces  +  X  StructureForces  (22a) 

{ms  +  Mf)vn  =  ^  Fluid  Forces + ^  StructureForces  (22b) 

where  ms  is  the  solid  mass  of  the  node,  mr  is  the  fluid 
mass  of  the  node,  im  and  yn  are  respectively  the  fluid 
and  solid  accelerations  of  the  node  normal  to  the  fluid- 
structure  interface.  The  tangential  fluid  and  solid 
accelerations  (ut,vt)  are  calculated  using  the  following 
equations: 

((1  -  s)ms  +  mf)iit  =  (1  -  s)^  StructureForces  +  ^  Fluid  Forces  (23a ) 

(ms  +  (1  -  s)m/)vt  =  ^  StructureForces  +  (1  -  s)^t  Fluid  Forces  (23  b) 

where  s  is  the  slip  factor.  A  no-slip  condition 
corresponds  to  a  slip  factor  of  zero.  The  slip  factor 
determines  how  much  of  the  fluid  and  structure  forces 
are  mutually  exchanged.  Equations  22  and  23  are 
written  for  all  fluid-solid  interface  nodes.  The  fluid  mesh 
must  move  and  deform  with  the  tank.  This  is  done  by 
modeling  the  fluid  mesh  using  very  light  and  compliant 
(3  orders  of  magnitude  less  than  the  tank)  solid  brick 
elements  (called  “mock”  mesh).  The  ALE  formulation  is 
used  to  account  for  the  fluid  mesh  deformation/motion. 

3.  TIRE  FRICTIONAL  CONTACT  MODEL 

A  tire  is  mounted  on  a  rigid  body  representing  the  wheel. 
Contact  is  detected  between  the  tire’s  surface  and  a 
polygonal  surface  representing  the  contact  body.  The 
tire’s  external  surface  is  discretized  in  the  circumference 
and  meridian  directions  into  a  grid  of  rectangles  with  a 
contact  point  defined  at  the  center  of  each  rectangle 
(Figure  3).  The  global  position  (  xgcp  )  and  velocity  ( xgcp  ) 
of  a  contact  point  is  obtained  using  the  rigid  body 
position  and  rotation  matrix: 

XGcPi  —  XbF-i  +  RBFy  XLcpj  (24) 

XGcpf  =  Xbfi  +  Ruff  (Wbfx  xlcP)j  (25) 

where  Xbf  and  Xbf  are  the  global  position  an  velocity 
vectors  of  the  wheel’s  frame,  Rbf  is  the  rotation  matrix 
of  the  wheel  relative  to  the  global  reference  frame,  Wbf 
is  the  wheel’s  angular  velocity  vector  relative  the  local 
wheel’s  frame,  and  xuP  is  the  position  of  the  contact 
point  relative  to  the  wheel’s  frame.  The  frictional  contact 
force  at  each  contact  point  (sum  of  the  normal  contact 
and  tangential  friction  forces)  is  transferred  as  a  force 
and  a  moment  to  the  center  of  the  wheel.  The  negative 
of  this  force  is  transferred  to  the  center  of  the  contact 


body.  Those  forces  and  moments  are  summed  and 
applied  as  external  forces  to  the  wheel  and  contact 
body. 


Figure  3  Tire  mounted  on  a  wheel  and  discritized  into  a  grid  of 
rectangular  elements  with  a  typical  contact  point  at  the  center  of  an 
element. 

3.1  PENALTY  NORMAL  CONTACT  MODEL 

The  penalty  technique  is  used  for  imposing  the  contact 
constraints.  In  this  technique,  a  normal  reaction  force 
(Fmrmai)  is  generated  when  a  contact  point  is  in  contact 
with  a  body  (i.e.  when  the  point  is  inside  the  body)  that  is 
given  by  [25,  26]: 

Fnormal  =  AkPd  +  A\  C”d .  d-0  (26) 

[, SpCpd  d  <  0 

where  A  is  the  area  of  the  rectangle  associated  with  the 
contact  point,  kP  and  cP  are  the  penalty  stiffness  and 
damping  coefficient  per  unit  area;  d  is  the  closest 
distance  between  the  node  and  the  contact  surface 
(Figure  4);  sP  is  a  separation  damping  factor  between  0 
and  1  which  determines  the  amount  of  sticking  between 
the  contact  node  and  the  contact  surface  at  the  node 
(leaving  the  body).  The  normal  contact  force  vector  is 
given  by: 

Fnt  —  /?,  F normal  (27) 

where  n  is  the  normal  to  the  surface.  The  total  force  on 
the  node  generated  due  to  the  frictional  contact  between 
the  point  and  surface  is  given  by: 

Fpo  inti  =  Ft-  +  Fni  (28) 


Contact 


Figure  4  Contact  surface  and  contact  node. 


3.2  ASPERITY  FRICTION  MODEL 

The  frictional  contact  force  (given  by  the  asperity  friction 
model  described  below)  transmitted  to  the  contact  body 
at  the  contact  point  is  Ft is  given  by: 

Ftf  =  Ftangent  ti  (29) 

An  asperity  friction  model  is  used  along  with  the  normal 
force  to  calculate  the  tangential  friction  force  (Ftangent)  [25]. 
The  basic  idea  of  the  model  is  that  friction  between  two 
rough  surfaces  in  contact  arises  due  to  the  interaction  of 
the  surface  asperities.  When  two  surfaces  are  in  static 
(stick)  contact,  the  surface  asperities  act  like  tangential 
springs.  When  a  tangential  force  is  applied,  the  springs 
elastically  deform  and  pull  the  surfaces  to  their  original 
position.  If  the  tangential  force  is  large  enough,  the 
surface  asperities  yield  (i.e.  the  springs  break)  allowing 
sliding  to  occur  between  the  two  surfaces.  The 
breakaway  force  is  proportional  to  the  normal  contact 
pressure.  In  addition,  when  the  two  surfaces  are  sliding 
past  each  other,  the  asperities  provide  resistance  to  the 
motion  that  is  a  function  of  the  sliding  velocity  and 
acceleration,  and  the  normal  contact  pressure.  Figure  5 
shows  a  schematic  diagram  of  the  asperity  friction 
model.  It  is  composed  of  a  simple  piece-wise  linear 
velocity-dependent  approximate  Coulomb  friction 
element  (that  only  includes  two  linear  segments)  in 
parallel  with  a  variable  anchor  point  spring. 


3.3  CONTACT  SEARCH 

Contact  between  each  contact  point  on  the  tire  and  the 
contact  surface(s)  (which  are  polygonal  surfaces)  is 
detected  using  a  binary  tree  contact  search  algorithm 
which  allows  fast  contact  search.  The  algorithm  works 
by  recursively  dividing  the  polygonal  surface  into  2 
blocks  of  polygons  and  then  finding  the  bounding  box  for 
each  block  of  polygons.  If  the  contact  point  is  inside  a 
bounding  box  then  the  two  sub-bounding  boxes  are 
checked  to  determine  if  the  point  is  inside  either  one. 
The  recursion  stops  when  the  bounding  box  contains  3 
polygons  or  less.  At  that  point  a  more  computationally 
intensive  contact  algorithm  between  a  point  and  a 
polygon  is  used  to  determine  the  depth  of  contact. 

4.  VOF  FREE-SURFACE  MODEL 

For  each  fluid  element  a  VOF  value  between  0  and  1  is 
defined,  where  0  corresponds  to  empty  elements  and  1 
corresponds  to  elements  completely  filled  with  fluid.  The 


elements’  VOF  values  are  updated  each  time-step  by 
moving  fluid  from  a  completely  or  partially  filled  “donor” 
element  to  an  empty  or  partially  filled  neighboring 
“acceptor”  element  using  the  following  model: 


Veo  =  Ve  VOF, 


Vna  =  Vn  (1-VOF,,) 


[A  t  S  A  n.u 


A.V  =  < 


Veo 

Vna 


(30) 

(31) 

Veo  >  AV  and  Vna  >  Veo 

Veo  <  AV 

(32) 

Vna  <  AV 

where  Ve  is  the  volume  of  the  element;  Vn  is  the  volume 
of  the  neighboring  element;  Veo  is  the  volume  of  the 
element  occupied  by  the  fluid;  Vna  is  the  volume  of  the 
neighboring  element  available  to  receive  fluid;  AV  is  the 
volume  flow  through  the  boundary  between  the  two 
elements  in  a  time  step;  At  is  the  solution  time  step;  S  is 
the  surface  area  between  the  two  elements;  A  is  a  value 
between  0  and  1  indicating  the  free-surface  aperture 
through  which  the  fluid  can  move  from  the  element  to 
the  neighboring  element;  n  is  a  unit  vector  normal  to  S; 
and  u  is  the  fluid  velocity  vector  at  the  surface  S.  If  AV  is 
less  than  0  then  the  element  is  an  acceptor  element  and 
the  VOF  values  are  not  updated  because  they  will  be 
updated  later  when  the  neighboring  element  is  set  to  be 
the  donor  element.  If  AV  is  greater  than  0  then  the  VOF 
values  are  updated  using  the  following  equations: 


VOF,  =  VOF,  -AV/Ve  (33) 

VOF»  =  VOF»  + AF  /Vn  (34) 


The  free-surface  apertures  A  at  the  element  interfaces 
are  used  to  limit  the  fluid  flow  based  on  the  location  of 
the  free  surface  inside  the  element.  A  is  calculated  as 
follows.  If  the  VOF  value  of  the  element  is  1  then  there  is 
no  free-surface  at  the  element,  therefore  A= 1.  For 
elements  with  a  VOF  value  less  than  1,  the  following 
steps  are  used  to  calculate  A: 


•  Calculate  the  normal  to  the  surface  by  looking  at  a 
stencil  of  neighboring  elements  around  the  element. 
This  is  done  using  the  following  equation: 

nei=VOFnkSnknnkl  /=  1 , 2,  3  (35) 

where  is  the  *  component  of  the  normal  to  the 
free-surface  at  the  element,  VOFb*  is  the  VOF  value 
for  neighboring  element  number  k,  Snk  is  the  area  of 
the  intersection  surface  between  the  element  and 
neighboring  element  k ,  and  nnki  is  the  component  i  of 
the  normal  to  the  surface  between  the  element  and 
neighboring  element  k.  He  is  then  normalized  into  a 
unit  vector.  Figure  6  shows  a  2D  4-node  quadrilateral 
and  the  free-surface  along  with  the  normal  ne  ■ 


•  Calculate  the  apertures  A  for  each  neighboring 
element  by  constructing  a  planar  surface  with 
normal  m  and  with  total  volume  equal  to  VOF,  Ve 
(see  Figure  6). 


Figure  6  Stencil  of  neighboring  elements  used  to  determine  the  free- 
surface  normal  ne  and  the  liquid  free-surface. 

5.  VALIDATION  STUDY 

The  model  is  validated  using  a  full-scale  army  heavy 
class  tactical  PLS  (Palletized-Load  System)  trailer 
carrying  a  potable  water  tank  module  (called  Hippo).  The 
trailer  was  placed  on  an  n-post  motion  base  simulator  in 
TARDEC’s  Simulation  Laboratory  (TSL)  (Figure  7).  The 
n-post  motion  simulator  consists  of  linear  hydraulic 
actuators  each  placed  under  one  of  the  trailer’s  tires. 
Each  actuator  has  one  degree  of  freedom  along  the 
vertical  direction  and  can  be  independently  commanded 
to  follow  a  certain  vertical  displacement  time-history. 
Thus,  the  actuators  can  be  controlled  in  such  a  way  as 
to  simulate  vertical  position  time-histories  of  the  wheels 
during  typical  road  maneuvers.  These  maneuvers  can 
include:  traversing  a  bumpy  terrain,  going  over 
symmetric  or  asymmetric  bumps,  turning  and  lane- 
change.  Note  that  the  motion  simulator  cannot  simulate 
the  inertial  centrifugal  and  Coriolis  forces  that  arise  due 
to  the  time-varying  motion  of  the  trailer  on  the  road. 


Figure  7  Experimental  setup  consisting  of  a  PLS  trailer  carrying  a 
“Hippo”  tank  module  mounted  an  n-post  motion  simulator  and 
instrumented  to  measure  the  tire  and  suspension  system  deflection 
time-histories. 

The  trailer  in  our  experimental  setup  has  3  axles.  We 
call  the  axles:  front  axle,  back  axle  1  and  back  axle  2 
(Figure  8).  The  trailer  also  has  6  wheels:  front  right,  front 
left,  back  1  right,  back  1  left,  back  2  right  and  back  2  left. 

Each  hydraulic  actuator  is  composed  of  a  cylinder, 
piston  and  a  dishpan  where  the  tire  rests.  The  front  and 
back  axle  1  actuators  had  trapezoidal  dishpans  for 


safety  reasons.  This  is  in  order  to  ensure  that  the  trailer 
does  not  slide  off  the  motion  simulator  (Figure  9).  The 
two  back  axle  2  actuators  had  flat  dishpans.  In  addition, 
a  wire  harness  was  attached  to  the  top  of  the  trailer  at  all 
time  during  the  experiments  for  safety  reasons.  The  wire 
harness  was  loose  so  as  not  to  interfere  with  the 
vibration  results  of  the  trailer. 


Figure  8  Bottom  view  of  the  trailer  showing  the  three  axles. 
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Figure  9  Side  view  of  the  trailer  (left)  showing  the  trapezoidal  dishpans 
in  yellow  on  the  font  axle  tires  and  back  axle  1  tires.  Back  axle  2  tires 
have  flat  dishpans. 


Displacement  data  is  measured  using  LVDTs.  All  LVDT 
signals  are  sampled  at  a  rate  of  256  samples/sec. 
Figure  10  shows  a  picture  of  the  trailer  along  with  the 
displacement  data  collected  during  the  experiments.  The 
experiment  displacement  measurements  are  compared 
with  the  results  predicted  using  the  computational 
model.  Displacement  is  measured  at  22  points.  They 
are: 

•  Linear  vertical  input  motion  of  the  6  actuators 
(labeled:  Front  left  actuator;  Front  right  actuator; 
Back  1  left  actuator;  Back  1  right  actuator;  Back  2 
left  actuator;  Back  2  right  actuator). 

•  Linear  vertical  motion  of  each  wheel’s  center  relative 
to  the  trailer’s  frame  (labeled:  Front  left  wheel;  Front 
right  wheel;  Back  1  left  wheel;  Back  1  right  wheel; 
Back  2  left  wheel;  Back  2  right  wheel). 

•  Linear  vertical  deflection  of  each  tire  obtained  by 
measuring  the  distance  between  the  center  of  the 
wheel  and  the  dishpan  (labeled:  Front  left  tire;  Front 
right  tire;  Back  1  left  tire;  Back  1  right  tire;  Back  2  left 
tire;  Back  2  right  tire). 

•  Longitudinal  displacements  of  the  trailer’s  frame 
relative  to  ground  measured  on  the  Left  and  right  of 
the  trailer’s  frame  (Labeled:  Left  long.;  Right  long.). 

•  Lateral  displacements  of  the  trailer’s  frame  relative 
to  ground  measured  near  the  front  and  back  2  axles 
(Labeled  Front  lateral;  Back  lateral). 


A  camera  mounted  on  the  ground  is  used  to  record  the 
motion  of  the  trailer  tank.  The  camera  is  set  to  capture 
30  frames/sec.  The  camera  view  is  shown  in  Figure  7. 


Figure  10  Displacement  data  collected  during  the  experiments. 


Table  2  Trailer  empty  tank  experiments  (26  experiments). 


Pitch 

Roll 

Stir 

Freq. 

(Hz) 

Amp. 

(mm) 

Freq. 

(Hz) 

Amp. 

(mm) 

Freq. 

(Hz) 

Amplitude 

(mm) 

Harmonic 

excitation 

0.7 

100 

0.5 

36 

0.5 

50 

0.7 

140 

0.5 

48 

0.5 

70 

1.0 

60 

1.0 

20 

0.6 

50 

1.0 

75 

1.0 

30 

0.6 

70 

1.5 

20 

1.5 

5 

1.0 

30 

1.5 

25 

1.5 

7 

1.0 

40 

2.0 

17 

2.0 

28 

Ramp 

excitation 

Amp.  (mm) 

Amplitude  (mm) 

Amplitude  (mm) 

63  (0.1  sec  ramp) 

45  (0.1  sec  ramp) 

155  (O.lsec  ramp  LF) 

110(0.2sec  ramp) 

50  (0.2sec  ramp) 

110  (O.lsec  ramp  LB) 

Table  3  Trailer  65%  filled-tank  experiments  (Total  of  26  experiments). 


Pitch 

Roll 

Stir 

Freq. 

(Hz) 

Amp. 

(mm) 

Freq. 

(Hz) 

Amp. 

(mm) 

Freq.  (Hz) 

Amp. 

(mm) 

Harmonic 

excitation 

0.7 

100 

0.5 

36 

0.5 

50 

0.7 

140 

0.5 

48 

0.5 

70 

1.0 

60 

1.0 

20 

0.6 

50 

1.0 

75 

1.0 

30 

0.6 

70 

1.5 

20 

1.5 

8 

1.0 

30 

1.5 

25 

1.5 

12 

1.0 

40 

2.0 

17 

2.0 

30 

Ramp 

excitation 

Amp.  (mm) 

Amp.  (mm) 

Amplitude  (mm) 

70  (O.lsec  ramp) 

50  (0.1  sec  ramp) 

180  (O.lsec  ramp  LF) 

120  (0.2sec  ramp) 

55  (0.2  sec  ramp) 

110  (O.lsec  ramp  LB) 

Experiments  were  carried  out  with  an  empty  tank  and  a 
65%-filled  tank.  The  time-histories  of  the  displacements 
at  the  22  points  described  above  are  collected  for 
various  input  motion  excitations.  Those  include  pitch,  roll 
and  stir.  Pitch  excitation  simulates  both  sides  of  the 
trailer  going  over  a  symmetric  bump/pothole.  Roll  motion 
simulates  the  trailer  turning  or  going  over  a 
bump/pothole  only  on  one  side  of  the  trailer.  The  stir 
motion  simulates  a  combination  of  pitch  and  roll  motions. 
For  each  motion  type  either  a  harmonic  excitation  or  a 
ramp  excitation  was  used.  For  the  harmonic  excitation 
the  frequency  and  amplitude  of  the  excitation  is  varied. 
The  harmonic  excitation  experiments  are  conducted  by 
slowly  ramping  up  the  amplitude  of  the  excitation  to  the 
full  amplitude.  For  the  ramp  excitation  just  the  amplitude 
of  the  excitation  is  varied.  Table  2  shows  the 
experiments  carried  out  with  the  empty  tank.  Table  3 
shows  the  experiments  carried  out  with  the  tank  65% 


filled  with  water.  For  the  harmonic  excitations  the 
specified  amplitude  is  the  peak-to-peak  amplitude.  For 
the  ramp  excitations,  the  amplitude  is  the  difference 
between  the  initial  value  and  the  final  value  after  the 
ramp.  The  two  actuators  under  the  left  two  back  tires  are 
always  moved  the  same  way.  Similarly,  the  two 
actuators  under  the  right  two  back  tires  are  always 
moved  the  same  way.  The  experiments  are  performed 
as  follows: 

•  In  the  pitch  experiments  the  left  and  right  actuators 
move  in  phase.  The  actuators  for  the  two  back  axles 
are  moved  180°  out-of-phase  with  the  front 
actuators. 

•  In  the  roll  experiments  left  and  right  actuators  move 
1 80°  out-of-phase.  The  two  back  axles  and  the  front 
actuators  move  in-phase. 

•  In  the  harmonic  stir  experiments  the  front  left,  front 
right,  back  right,  and  back  left  actuators  are  90°  out- 
of-phase. 

•  In  the  ramp  stir  experiments  either  the  left  front  (LF) 
or  the  two  left  back  (LB)  actuators  are  moved. 


Figure  11  Trailer  model. 


Figure  12  Front  leaf-spring  (top)  back  leaf-spring  (bottom)  modeled 
using  brick  elements.  Connection  points  are  shown  as  green  spheres. 

The  modeling  techniques  presented  in  this  paper  were 
implemented  in  the  DIS  [32]  commercial  finite  element 
code.  The  DIS  code  was  used  to  numerically  simulate 


the  experiments  in  Tables  2-3.  The  simulation  results 
are  then  compared  to  the  experiment  results.  A  DIS 
finite  element  model  of  the  trailer  and  the  n-post  motion 
simulator  was  constructed.  The  model  has  33548  nodes 
and  consists  of  the  following  components: 

•  A  rigid  grounded  base. 

•  Rigid  bodies  representing  the  chassis,  3  axles,  and 
suspension  system  elements  (Figure  11). 

•  6  linear  actuators.  The  actuator  cylinders  are  part  of 
the  grounded  base  and  the  actuator  pistons  and 
dishpans  are  modeled  as  rigid  bodies  (Figure  10). 
Two  parallel  cylindrical  joints  are  used  at  each 
actuator  to  model  a  prismatic  joint  that  allows  the 
actuator  to  move  only  along  the  vertical  direction. 

•  2  front  and  2  back  leaf-springs  modeled  using  brick 
elements  (Figure  12). 

•  6  wheels  modeled  using  rigid  bodies.  A  tire  is 
mounted  on  each  wheel.  The  tire  is  in  contact  with 
the  dishpan  of  the  corresponding  actuator  piston 
(Figure  9). 

•  An  oval  tank.  The  tank  is  modeled  as  a  rigid  body. 
The  tank  is  discretized  using  33516  hexahedral  fluid 
elements  (Figure  13).  The  tank  has  a  cross-section 
baffle  near  it’s  center.  The  baffle  has  a  big  round 
opening  at  the  center  and  small  openings  near  the 
bottom  and  top  of  the  tank  to  equalize  the  liquid 
level. 

•  Spherical  joints  are  used  to  model  the  suspension 
system  connections  and  the  wheels’  connections  to 
the  axles. 


Figure  13  Tank  geometric  model  (top)  finite  element  model  (bottom). 


Water  (p  =  1000  Kg/m3,  //=  0.001  Kg/m. sec)  is  modeled 
as  incompressible  using  the  artificial  compressibility 
technique  with  an  artificial  sound  speed  factor  of  0.1  (i.e. 
the  artificial  sound  speed  in  the  water  is  taken  as  1483 
m/sec  x  0.1  =  148.3  m/sec).  Due  to  the  use  of  large 
elements  near  the  solid  surface,  full  slip  boundary 
condition  at  the  wall  is  used  (s  =  1  in  Equation  23).  Thus, 
the  viscous  wall  friction  effects  are  assumed  to  be 
negligible.  Gravity  is  modeled  with  the  gravitational 
acceleration  taken  to  be  9.8  m/sec2  in  the  vertical 
direction.  The  explicit  time  step  was  1 .53x1 0"5  sec. 

After  performing  the  simulations,  we  determined  that  the 
trapezoidal  dishpans  that  were  used  for  the  front  and 
back  axle  1  tires  in  order  to  ensure  that  the  trailer  does 
not  slide  off  the  motion  simulator,  affected  the 
displacement  response  at  those  dishpans.  This  is  due  to 
the  fact  that  the  sides  of  the  tires  were  in  contact  with 
the  sides  of  the  trapezoidal  dishpans.  This  produced 
friction  forces  which  made  the  tire  stick  to  the  dishpan 
whenever  the  dishpan  is  trying  to  pull  away  from  the  tire. 
The  1-node  tire  model  is  not  adequate  for  modeling  this 
because  it  does  not  account  for  the  tire  deformation  with 
respect  to  the  wheel.  This  is  due  to  the  fact  that  the  finite 
element  node  of  the  wheel  is  also  used  for  the  tire 
contact  and  deformation  calculations.  In  order  to  account 
for  the  deformation  of  the  tire  with  respect  to  the  wheel  a 
detailed  finite  element  model  of  the  tire  is  needed  similar 
to  the  model  developed  in  Ref.  [27],  For  practical 
applications,  the  side  of  the  tire  is  not  in  contact  with  the 
terrain,  therefore,  the  1-node  tire  model  can  model  with 
adequate  accuracy  the  tire’s  response.  The  flat  dishpans 
used  for  the  back  axle  2  tires  do  not  have  this  problem 
and  therefore  the  1-node  tire  model  is  adequate.  In  the 
present  paper,  all  the  runs  were  performed  using  the  1- 
node  tire  model.  Therefore  we  will  present  the  results  for 
the  back  axle  2.  In  a  future  paper,  we  will  present  the 
results  with  a  detailed  finite  element  tire  model  for  the 
trailer’s  tires. 

The  graphs  in  Figures  14-29  show  a  comparison  of  the 
DIS  simulation  and  experiment  displacements  at  the 
back  axle  2  of  the  trailer  (back  2  left  wheel;  back  2  right 
wheel;  back  2  left  tire;  back  2  right  tire)  for  eight  empty 
tank  and  eight  65%-filled  tank  experiments.  The  average 
difference  in  magnitude  and  frequency  between  the 
experiment  and  simulation  is  about  15-20%  in  the 
graphs  shown  in  Figures  14-29.  The  average  difference 
in  magnitude  and  frequency  between  the  experiment 
and  simulation  for  the  front  and  back  axle  2 
displacements  was  about  30%  on  average  due  to  the 
presence  of  the  trapezoidal  dishpans  as  mentioned 
above.  There  was  no  significant  difference  in  the 
magnitude  of  the  experiment  to  simulation  error  between 
the  empty  tank  runs  and  the  65%-filled  tank  runs.  In  the 
opinion  of  the  authors,  the  main  sources  of  error 
between  the  experiments  and  simulation  in  order  of 
importance  are: 

•  The  trapezoidal  dishpans  at  the  front  and  back  axle 
1  tires,  affected  the  back  axle  2  response. 


Figure  14  Comparison  between  the  experiment  and  DIS  simulation  for 


Figure  16  Comparison  between  the  experiment  and  DIS  simulation  for 
an  empty  tank  with  a  pitch  2.0  Hz,  17  mm  harmonic  excitation. 


Figure  15  Comparison  between  the  experiment  and  DIS  simulation  for 


Figure  17  Comparison  between  the  experiment  and  DIS  simulation  for 
a  65%  filled  tank  with  a  pitch  2.0  Hz,  17  mm  harmonic  excitation. 


Figure  18  Comparison  between  the  experiment  and  DIS  simulation  for 


Figure  20  Comparison  between  the  experiment  and  DIS  simulation  for 
an  empty  tank  with  a  roll  1 .0  Hz,  30  mm  harmonic  excitation. 


Figure  19  Comparison  between  the  experiment  and  DIS  simulation  for 


Figure  21  Comparison  between  the  experiment  and  DIS  simulation  for 
a  65%  filled  tank  with  a  roll  1.0  Hz,  30  mm  harmonic  excitation. 


Figure  22  Comparison  between  the  experiment  and  DIS  simulation  for 


Figure  24  Comparison  between  the  experiment  and  DIS  simulation  for 
an  empty  tank  with  a  stir  1 .0  Hz,  40  mm  harmonic  excitation. 


Figure  23  Comparison  between  the  experiment  and  DIS  simulation  for 


Figure  25  Comparison  between  the  experiment  and  DIS  simulation  for 
a  65%-filled  tank  with  a  stir  1.0  Hz,  40  mm  harmonic  excitation. 


Figure  26  Comparison  between  the  experiment  and  DIS  simulation  for 
an  empty  tank  with  a  pitch  0.1  sec.  63  mm  ramp  excitation. 
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Figure  27  Comparison  between  the  experiment  and  DIS  simulation  for 


Figure  28  Comparison  between  the  experiment  and  DIS  simulation  for 
an  empty  tank  with  a  roll  0.2  sec.  50  mm  ramp  excitation. 


Figure  29  Comparison  between  the  experiment  and  DIS  simulation  for 
a  65%-filled  tank  with  a  roll  0.2  sec.  55  mm  ramp  excitation. 
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Figure  30  Snapshots  of  the  trailer  and  liquid  sloshing  for  the  pitch  0.7  Hz,  140  mm  harmonic  excitation  experiment. 
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Figure  31  Snapshots  of  the  trailer  and  liquid  sloshing  for  the  roll  0.5  Hz,  48  mm  harmonic  excitation  experiment. 


•  The  tank  module  (Hippo)  was  not  securely  mounted 
on  the  trailer.  There  was  a  clearance  in  the 
connection  between  the  Hippo  and  the  trailer,  which 
based  on  visual  observation  and  the  sound 
produced  is  probably  +/-  20  mm.  This  means  that 
the  Hippo  was  in  frictional  contact  with  the  trailer. 
For  the  low  frequency  runs  the  friction  was  enough 
to  lock  the  tank  on  the  trailer.  But  for  the  higher 
frequency  runs  and  for  the  ramp  runs  the  Hippo  was 
sliding  on  the  trailer.  In  the  simulation  model,  it  was 
assumed  that  the  connection  between  the  Hippo  and 
the  trailer  was  rigid. 

•  Tire  damping  as  a  function  of  tire  deflection  was 
estimated.  However,  this  damping  may  be  non¬ 
linear. 

•  Clearances  in  the  trailer  joints  were  not  modeled 
since  they  were  not  measured. 

•  Non-linear  behavior  of  the  suspension  leaf-springs. 
Friction  and  damping  were  estimated  and  added  to 
the  leaf-springs  using  truss  elements.  It  is  very  hard 
to  accurately  characterize  the  non-linear  stiffness, 
damping  and  friction  of  the  leaf-springs. 

•  The  security  harness  on  the  tank  module  on  top  of 
the  PLS  trailer  contributed  to  the  difference  between 
the  simulation  and  experiments  especially  for  the 
higher  frequency  excitations. 

Figures  30  and  31  show  snapshots  of  the  simulation  of 
the  trailer  motion  and  liquid  sloshing  in  the  tank  for  two 
typical  experiments. 

6.  CONCLUDING  REMARKS 

A  finite  element  model  for  predicting  the  fully  coupled 
dynamic  response  of  flexible  multibody  systems  and 
liquid  sloshing  for  tanker  trucks  was  presented.  The 
model  has  the  following  characteristics: 

•  Parallel  explicit  time-integration  solver. 

•  Library  of  accurate  large  rotation  finite  elements 
including:  truss,  beam,  shell  and  solid  elements.  The 
elements  only  use  Cartesian  coordinates  as  DOFs. 

•  The  fluid  mesh  is  modeled  using  a  very  light  and 
compliant  solid  mesh  which  allows  the  fluid  mesh  to 
move/deform  along  with  the  tank  using  the  Arbitrary 
Lagrangian-Eulerian  formulation. 

•  Acceptor-donor  VOF  algorithm  for  modeling  the 
fluid's  free-surface. 

•  The  motion  of  the  solid  and  fluid  is  referred  to  a 
global  inertial  Cartesian  reference  frame. 

•  A  total  Lagrangian  deformation  description  is  used 
for  the  solid  elements. 

•  The  penalty  technique  is  used  to  model  the  joints. 

•  1 -Node  tire  model. 

A  validation  study  of  the  finite  element  model  was 
carried  out  using  a  tanker-trailed  mounted  on  an  n-post 


motion  simulator.  The  study  shows  that  the  model  can 
predict  reasonably  well,  within  15-20%  difference  on 
average,  the  response  of  the  trailer.  There  was  no 
significant  difference  in  error  magnitude  between  the 
simulation  and  experiment  for  the  empty  tank  and  65% 
filled  tank  runs.  This  shows  that  the  difference  between 
the  simulation  and  experiment  is  not  due  to  the  fluid- 
structure  interaction  modeling  and  is  mostly  due  to  the 
trailer  multibody  model  error  sources  identified  in 
Section  5. 
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APPENDIX:  EXPLICIT  SOLUTION  PROCEDURE 

The  solution  fields  for  modeling  the  solid,  fluid  and  liquid 

free-surface  are  defined  at  the  model  nodes.  These  are: 

•  Solid  translational  positions. 

•  Solid  translational  velocities. 

•  Solid  translational  accelerations. 

•  Solid  rotation  matrices. 

•  Solid  rotational  (angular)  velocities. 

•  Solid  rotational  (angular)  accelerations. 

•  Fluid  velocities. 

•  Fluid  accelerations. 

•  Fluid  pressure. 

•  Fluid  pressure  rate. 

•  Volume-of-fluid. 


•  Eddy  kinetic  energy 

•  Eddy  kinetic  energy  rate. 

The  explicit  time  integration  solution  procedure  for 
modeling  the  coupled  response  of  the  solid  (multibody 
system),  fluid,  and  liquid  free-surface  (using  the  VOF 
formulation)  predicts  the  time  evolution  of  the  above 
response  quantities.  The  procedure  is  implemented  in 
the  DIS  [32]  (Dynamic  Interactions  Simulator) 
commercial  software  code  and  is  outlined  below: 

1)  Prepare  the  run: 

a.  Set  the  initial  conditions  for  the  solution  fields 
identified  above. 

b.  Create  a  list  of  all  the  finite  elements  (including 
both  solid  and  fluid  elements). 

c.  Create  a  list  of  elements  that  will  run  on  each 
processor.  This  is  done  using  an  algorithm  which 
tries  to  make  the  computational  cost  on  each 
processor  equal. 

d.  Create  a  list  of  all  the  constraints  (including  both 
solid  and  fluid  constraints). 

e.  Calculate  the  solid  masses  for  each  finite  element 
node  by  looping  through  the  list  of  finite  elements. 
Note  that  the  solid  masses  are  fixed  in  time. 

f.  For  each  node  create  a  list  of  corner  and  edge 
nodes  that  are  connected  to  it  using  fluid  volume 
elements. 

g.  VOF  preparations: 

i.  Find  a  list  of  the  volume  fluid  elements. 

ii.  Create  a  list  of  fluid  volume  elements  that  will 
run  on  each  processor.  This  is  done  using  an 
algorithm  which  tries  to  make  the  computational 
cost  on  each  processor  equal. 

iii.  For  each  element  find  all  neighboring  elements. 

iv.  For  each  element  find  the  element  VOF  using 
the  nodal  VOF. 

v.  Re-interpolate  the  elements’  VOF  to  nodal  VOF. 

h.  Loop  over  all  the  elements  and  find  the  minimum 
time  step  for  the  explicit  solution  procedure. 

i.  Loop  over  all  the  elements  and  create  a  list  of  wall 
nodes.  For  each  wall  node  find  the  list  of  fluid 
boundary  elements. 

2)  Loop  over  the  solution  time  and  increment  the  time 
by  At  each  step  while  doing  the  following: 

a.  Set  the  nodal  values  at  the  last  time  step  to  be 
equal  to  the  current  nodal  values  for  all  solution 
fields. 

b.  Do  2  iterations  (a  predictor  iteration  and  a 
corrector  iteration)  of  the  following: 

i.  Initialize  the  nodal  fluxes  to  zero.  Those  include: 
solid  forces,  solid  moments,  fluid  forces, 
boundary  fluid  forces,  and  pressure  fluxes.  In 
addition,  the  lumped  nodal  fluid  volume  and  fluid 
mass  vectors  are  also  initialized  to  zero. 

ii.  Calculate  the  nodal  solid  and  fluid  fluxes  and  the 
lumped  fluid  volume/mass  vectors  by  looping 
through  all  the  elements  while  calculating  and 
assembling  the  element  nodal  fluxes  and 
vectors.  This  is  the  most  computationally 
intensive  step.  This  step  is  done  in  parallel  by 


running  each  list  of  elements  identified  in  step 
1  .c  on  one  processor. 

iii.  Find  the  nodal  values  at  the  current  time  step 
using  the  semi-discrete  equations  of  motion  and 
the  trapezoidal  time  integration  rule  (Equations 
1-5  and  19-21). 

iv.  Execute  the  solid  and  fluid  constraints.  The 
constraints  prescribe  the  nodal  values. 

v.  Apply  fluid-structure  interface  boundary 
conditions  for  all  wall  nodes  found  in  Step  l.i 
(see  Equations  22,  23).  This  is  done  by  doing 
the  following  for  each  wall  node: 

1.  Find  the  normal  to  the  surface  at  the  wall 
node. 

2.  Normalize  the  surface  normal. 

3.  Find  the  solid,  fluid  bulk  and  fluid  boundary 
forces  in  the  directions  normal  and  tangent  to 
the  surface. 

4.  Find  the  normal  and  tangential  solid  and  fluid 
accelerations  using  the  trapezoidal 
integration  rule  and  the  wall  slip  percentage. 

vi.  Set  the  pressure  boundary  conditions  at  the  free 
surface. 

vii.  Update  the  VOF  field: 

1 .  For  each  fluid  element  calculate  the  element 
volume.  This  step  is  done  in  parallel  using 
the  list  of  fluid  elements  for  each  processor 
found  in  Step  l.g.ii. 

2.  For  each  fluid  element  find  the  apertures 
through  which  the  fluid  convects  to  each 
neighboring  element.  This  step  is  done  in 
parallel  using  the  list  of  fluid  elements  for 
each  processor  found  in  Step  l.g.ii. 

3.  For  each  fluid  element  use  the  apertures,  the 

element  volume,  the  element  current  VOF 
value,  and  the  element  nodal  velocities  to 
update  the  VOF  value  of  all  neighboring 
elements  by  finding  the  volume  of  fluid  that 
left  the  element  during  that  time  step  using 
Equations  30-35.  This  step  is  done  in  parallel 
using  the  list  of  fluid  elements  for  each 
processor  found  in  Step  l.g.ii.  Note  that  this 
step  depends  on  the  order  of  the  elements  in 
the  list  of  elements.  However,  since  the 
updates  of  the  VOF  field  between  solution 
time  steps  are  small,  therefore  this 

dependence  is  generally  very  small.  In  order 
to  assure  minimum  dependence  on  the 
elements’  order,  at  a  time  step  the  elements 
are  updated  from  first  to  last,  then  at  the  next 
time  step  they  are  updated  from  last  to  first. 

viii.  Average  the  fluid  pressure  (This  step  eliminates 
the  pressure  checker-boarding  effect  and  allows 
use  of  equal  order  interpolation  for  both 
pressure  and  velocity). 

ix.  Go  to  the  beginning  of  step  2. 

An  advantage  of  explicit  solution  procedures  is  that  they 
are  “embarrassingly”  parallel.  The  above  procedure 
achieves  near  linear  speed-up  with  the  number  of 
processors. 


